add single precision versions of linear/gamma conversions

diff --git a/babl/base/pow-24.c b/babl/base/pow-24.c
index a3bb36e32da999296136a5b6ddcb6d7dbb949df8..03fca4346e585f877593940dc05ace29d4b12747 100644 (file)
--- a/babl/base/pow-24.c
+++ b/babl/base/pow-24.c
@@ -71,3 +71,57 @@ babl_pow_1_24 (double x)
     y = (7./6.) * y - z * ((y*y)*(y*y)*(y*y*y));
   return x*y;
 }
+
+//////////////////////////////////////////////
+/* a^b = exp(b*log(a))
+ *
+ * Extracting the exponent from a float gives us an approximate log.
+ * Or better yet, reinterpret the bitpattern of the whole float as an int.
+ *
+ * However, the output values of 12throot vary by less than a factor of 2
+ * over the domain we care about, so we only get log() that way, not exp().
+ *
+ * Approximate exp() with a low-degree polynomial; not exactly equal to the
+ * Taylor series since we're minimizing maximum error over a certain finite
+ * domain. It's not worthwhile to use lots of terms, since Newton's method
+ * has a better convergence rate once you get reasonably close to the answer.
+ */
+static inline float
+init_newtonf (float x, float exponent, float c0, float c1, float c2)
+{
+    int iexp;
+    float y = frexpf(x, &iexp);
+    y = 2*y+(iexp-2);
+    c1 *= M_LN2*exponent;
+    c2 *= M_LN2*M_LN2*exponent*exponent;
+    return y = c0 + c1*y + c2*y*y;
+}
+
+/* Returns x^2.4 == (x*(x^(-1/5)))^3, using Newton's method for x^(-1/5).
+ */
+float
+babl_pow_24f (float x)
+{
+  float y = init_newtonf (x, -1.f/5, 0.9953189663f, 0.9594345146f, 0.6742970332f);
+  int i;
+  for (i = 0; i < 3; i++)
+    y = (1.f+1.f/5)*y - ((1./5)*x*(y*y))*((y*y)*(y*y));
+  x *= y;
+  return x*x*x;
+}
+
+/* Returns x^(1/2.4) == x*((x^(-1/6))^(1/2))^7, using Newton's method for x^(-1/6).
+ */
+float
+babl_pow_1_24f (float x)
+{
+  float y = init_newtonf (x, -1.f/12, 0.9976800269f, 0.9885126933f, 0.5908575383f);
+  int i;
+  float z;
+  x = sqrtf (x);
+  /* newton's method for x^(-1/6) */
+  z = (1.f/6.f) * x;
+  for (i = 0; i < 3; i++)
+    y = (7.f/6.f) * y - z * ((y*y)*(y*y)*(y*y*y));
+  return x*y;
+}

diff --git a/babl/base/pow-24.h b/babl/base/pow-24.h
index e82229753796473056c718a7a8f267e53ede68b1..2fb338fe8e7b492a6021cbc55f1717d5f69a4460 100644 (file)
--- a/babl/base/pow-24.h
+++ b/babl/base/pow-24.h
@@ -21,5 +21,7 @@
 
 extern double babl_pow_1_24 (double x);
 extern double babl_pow_24 (double x);
+extern float  babl_pow_1_24f (float x);
+extern float  babl_pow_24f (float x);
 
 #endif

diff --git a/babl/base/util.h b/babl/base/util.h
index cfb17b0878e3c82822c448e44845b522262f29ae..34725230274738d7b4072792e2aa8ede2f5d1b07 100644 (file)
--- a/babl/base/util.h
+++ b/babl/base/util.h
@@ -78,6 +78,14 @@ babl_linear_to_gamma_2_2 (double value)
     return 1.055 * babl_pow_1_24 (value) - 0.055;
   return 12.92 * value;
 }
+static inline float
+babl_linear_to_gamma_2_2f (float value)
+{
+  if (value > 0.003130804954f)
+    return 1.055f * babl_pow_1_24f (value) - 0.055f;
+  return 12.92f * value;
+}
+
 
 static inline double
 babl_gamma_2_2_to_linear (double value)
@@ -86,6 +94,13 @@ babl_gamma_2_2_to_linear (double value)
     return babl_pow_24 ((value + 0.055) / 1.055);
   return value / 12.92;
 }
+static inline float
+babl_gamma_2_2_to_linearf (float value)
+{
+  if (value > 0.04045f)
+    return babl_pow_24f ((value + 0.055f) / 1.055f);
+  return value / 12.92f;
+}
 
 #else
   #define linear_to_gamma_2_2(value) (pow((value), (1.0F/2.2F)))